# ref to: http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/116-mfa-multiple-factor-analysis-in-r-essentials/

# load libs
library("FactoMineR")
library("factoextra")
library("dplyr")
library("gridExtra")
library("reshape2")
library("corrplot")
library("rmarkdown")
library("rgl")
library("NbClust")

config.rodionov.L <- 7 # rodionov regime-shift minimum length (7 = regime shift can happens not ofter then 1 time in 7 years)
config.rodionov.alpha <- 0.05 # confidence level, alpha = 0.05 equals p = 0.95

source("rodionov.R")

# load data
data <- read.csv("input/azov-dataset.csv")
# subset starts with 1972
data <- data[which(data$year >= "1972"),]


# helper functions
plot.anomaly <- function(x, y, title) {
  r <- barplot(y, col = "black", space = 0.7, 
               axis.lty = "solid", 
               names.arg = x, axisnames = TRUE, 
               ylab = "Anomalies", 
               main = title)
  box(lty = "solid", col = "black")
  # creates a horizontal line at y = 0
  abline(h = 0)
}

# step 1. Look up for anomaly diagnostics
data.anomaly <- data 
for (i in 2:ncol(data)) {
  mean <- mean(data.anomaly[,i], na.omit =T)
  data.anomaly[,i] <- data.anomaly[,i] - mean(data.anomaly[,i])
}

# save to png plots: environment summary (o2, sst, ppn, pcb, phosphate)
png("output/anomaly/env.png", width=1900, height=1600, res=200, type="cairo")
# plot 3 columns in 1 row
par(mfrow = c(2,2))
# plot anomalies
plot.anomaly(data.anomaly$year, data.anomaly$don_flow, "Don river flow")
plot.anomaly(data.anomaly$year, data.anomaly$salinity_taganrogbay, "Salinity Taganrog Bay")
plot.anomaly(data.anomaly$year, data.anomaly$salinity_azovsea, "Salinity Azov Sea")
plot.anomaly(data.anomaly$year, data.anomaly$temp_summer, "Summer water temp")
dev.off()

# save to png plots: food chain, consumers 1st rank
png("output/anomaly/food-consumers-1strank.png", width=1900, height=800, res=200, type="cairo")

par(mfrow = c(1,3))

plot.anomaly(data.anomaly$year, data.anomaly$zooplankton, "Zooplankton")
plot.anomaly(data.anomaly$year, data.anomaly$zoobentos, "Zoobentos")
plot.anomaly(data.anomaly$year, data.anomaly$mnemiopsis, "Mnemiopsis")

dev.off()

# save to png plots: pelagic consumers
png("output/anomaly/pelagic-consumers-2ndrank.png", width=1900, height=800, res=200, type="cairo")

par(mfrow = c(1,2))

plot.anomaly(data.anomaly$year, data.anomaly$clupionella_biomass, "Clupionella biomass")
plot.anomaly(data.anomaly$year, data.anomaly$anchovy_biomass, "Anchovy biomass")

dev.off()

# save to png plots: bottom consumers & anadromous fishes
png("output/anomaly/bottom-consumers-2ndrank.png", width=1900, height=1600, res=200, type="cairo")

par(mfrow = c(2,2))

plot.anomaly(data.anomaly$year, data.anomaly$gobidae_biomass, "Gobiidae biomass")
plot.anomaly(data.anomaly$year, data.anomaly$sander_biomass, "Sander biomass")
plot.anomaly(data.anomaly$year, data.anomaly$sturgeon_biomass, "Sturgeon biomass")
plot.anomaly(data.anomaly$year, data.anomaly$rutilus_biomass, "Rutilus biomass")

dev.off()

# save to png plots: human activities, anthropogenic impact
png("output/anomaly/human-impact.png", width=1900, height=1600, res=200, type="cairo")

par(mfrow = c(3,2))

plot.anomaly(data.anomaly$year, data.anomaly$sander_catches, "Sander catches")
plot.anomaly(data.anomaly$year, data.anomaly$gobidae_catches, "Gobiidae catches")
plot.anomaly(data.anomaly$year, data.anomaly$rutilus_catch, "Rutilus catches")
plot.anomaly(data.anomaly$year, data.anomaly$clupionella_catches, "Clupionella catches")
plot.anomaly(data.anomaly$year, data.anomaly$clupionella_catches, "Anchovy catches")


dev.off()

# save to png plots: human activities, anthropogenic impact
png("output/anomaly/liza-invade.png", width=1900, height=800, res=200, type="cairo")

par(mfrow = c(1,2))

plot.anomaly(data.anomaly$year, data.anomaly$liza_biomass, "Liza biomass")
plot.anomaly(data.anomaly$year, data.anomaly$liza_catches, "Liza catches")

dev.off()


# fit MFA model
data.clear <- data[,2:ncol(data)]
rownames(data.clear) <- data$year
# model info:
# --- groups: vector of factors, that assume related factors in 1 group.
# for this example, first 3 columns = sprat, next 5 columns = env, last 2 columns = plankton
# --- type: group variable type:
# s - quality variable, numerical with standartization
# c - quality variable, numerical NO standartization
# n - categorical variable
# f - frequency variable, 0 .. 1
fit <- MFA(data.clear, 
           group = c(4, 3, 2, 4, 5, 2), 
           type = c("s", "s", "s", "s", "s", "s"), 
           name.group = c("env", "food", "pelagic_bio", "bentic_bio", "anthrop.", "acclimat."),
)

# get summary individual impact
fit.impact <- get_mfa_var(fit, "quanti.var")
fit.score <- get_mfa_ind(fit)


# 1 - scree explained variance by dimensions
png("output/mfa/scree-plot.png", width=1900, height=1600, res=200, type="cairo")
fviz_screeplot(fit)
dev.off()

# 2 - factor & group & year contribution
png("output/mfa/factor-contrib-group.png", width=1900, height=1600, res=200, type="cairo")
print({grid.arrange(fviz_contrib(fit, choice = "group", axes = 1), 
             fviz_contrib(fit, choice = "group", axes = 2), 
             fviz_contrib(fit, choice = "group", axes = 3),
             nrow=3)
})
dev.off()

png("output/mfa/factor-contrib-var.png", width=1900, height=1600, res=200, type="cairo")
print({
  grid.arrange(fviz_contrib(fit, choice = "quanti.var", axes = 1), 
               fviz_contrib(fit, choice = "quanti.var", axes = 2), 
               fviz_contrib(fit, choice = "quanti.var", axes = 3),
               nrow=3)
})
dev.off()

png("output/mfa/factor-contrib-year.png", width=1900, height=1600, res=200, type="cairo")
print({
  grid.arrange(fviz_contrib(fit, choice = "ind", axes = 1), 
               fviz_contrib(fit, choice = "ind", axes = 2), 
               fviz_contrib(fit, choice = "ind", axes = 3),
               nrow=3)
})
dev.off()


# 3 - biplot main dim 1 vs dim2
png("output/mfa/biplot-dim1-vs-dim2.png", width=1900, height=1600, res=200, type="cairo")
print(fviz_mfa_var(fit))
dev.off()

# 4 - individual biplot for years
png("output/mfa/years-dim1-vs-dim2.png", width=1900, height=1600, res=200, type="cairo")
fviz_mfa_ind(fit, col.ind = "cos2", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))
dev.off()

# 5 - individuals with year trajectory of DIM1-2-3
png("output/mfa/biplot-years-traj.png", width=1900, height=1600, res=200, type="cairo")
fviz_mfa_ind(fit, partial = c("2022"))
dev.off()


# 6 - regime shift indication PC1-3
fit.score3 <- as.data.frame(fit.score$coord[,1:3])
years <- as.numeric(rownames(fit.score$coord))
fit.score3$year <- years
m <- melt(fit.score3, id="year")

png("output/mfa/preliminary-regime-shift.png", width=1900, height=1600, res=200, type="cairo")
print({ggplot(data=m, aes(x=year, y=value, colour=variable, size=variable, alpha=variable)) +
    geom_line(aes(linetype=variable)) + 
    geom_point(size = 2) + 
    scale_color_manual(values=c("red", "blue", "darkgreen"))+
    scale_size_manual(values=c(1.1, 1.0, 0.9),labels=c("Dim.1", "Dim.2", "Dim.3")) + 
    scale_alpha_manual(values=c(1,1,1)) + 
    xlab("Year") + 
    ylab("DIM variance") +
    geom_hline(yintercept = 0)})

dev.off()

# 8 - hierarhic clustering
hc <- hclust(dist(cbind(fit.score$coord[,1], fit.score$coord[,2], fit.score$coord[,3])), method = "ward.D2")
png("output/mfa/cluster-dendrogram.png", width=1900, height = 1600, res = 200, type="cairo")
plot(hc, main = "MFA cluster dendrogram by DIM variations impact", xlab = "Obs. year", ylab = "Cluster height")
rect.hclust(hc, k=6, border="gray")
dev.off()

### Perform K-means clustering procedure
# find optimum number of clusters
png("output/mfa/kmeans-cluster-silhouette.png", width=1900, height=1600, res=200, type="cairo")
print({fviz_nbclust(fit.score$coord, kmeans, method = "silhouette")})
dev.off()

# 6 clusters is optim
nb <- NbClust(fit.score$coord, distance = "euclidean", min.nc = 2, max.nc = 8, method = "kmeans")
clust.numb <- length(unique(nb$Best.partition))

kc12 <- kmeans(fit.score$coord[,c(1,2)], clust.numb)
kc13 <- kmeans(fit.score$coord[,c(1,3)], clust.numb)
kc23 <- kmeans(fit.score$coord[,c(2,3)], clust.numb)

png("output/mfa/kmeans-pc-1-2.png", width=1900, height=1600, res=200, type="cairo")
print({fviz_cluster(kc12, data = fit.score$coord[,c(1,2)], ggtheme = theme_minimal(), labelsize=16) +
    geom_hline(yintercept = 0, linetype = "dashed") + 
    geom_vline(xintercept = 0, linetype = "dashed")})
dev.off()

png("output/mfa/kmeans-pc-1-3.png", width=1900, height=1600, res=200, type="cairo")
fviz_cluster(kc13, data = fit.score$coord[,c(1,3)], ggtheme = theme_minimal()) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  geom_vline(xintercept = 0, linetype = "dashed")
dev.off()

png("output/mfa/kmeans-pc-2-3.png", width=1900, height=1600, res=200, type="cairo")
fviz_cluster(kc23, data = fit.score$coord[,c(2,3)], ggtheme = theme_minimal()) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  geom_vline(xintercept = 0, linetype = "dashed")
dev.off()


# 9 - traffic-light plot 

# get data dimensions
d <- dim(data)
data.names <- names(data)[-1]
# get dataset without year column
data.noyear <- data[,-1]

# prepare quantile statistics matrix
quant <- matrix(nrow = 5, ncol = (ncol(data.noyear)))
# fill matrix with quantile values for levels: 0.2, 0.4, 0.6, 0.8, 1
for (i in 1:ncol(data.noyear)) {
  quant[,i] <- quantile(data.noyear[,i], probs = c(0.2, 0.4, 0.6, 0.8, 1), na.rm = T)
}

# create quantile-relation matrix categorization
# focal point: create factor matrix with determination 
# to quantile intervals 
# ex, 1 means that value is in 0 ... 0.2 quantile of full factor row
quant.matrix <- data.noyear
for (i in 1:(nrow(quant.matrix))) {
  for (j in 1:(ncol(quant.matrix))) {
    cel <- quant.matrix[i, j]
    if (is.na(cel))
      quant.matrix[i,j] <- TRUE
    else if (cel < quant[1, j]) # cel in 1st quantile value between 0 ... 0.2 in col
      quant.matrix[i, j] <- 1 # set 1st quantile to output
    else if (cel >= quant[1, j] && cel < quant[2, j]) 
      quant.matrix[i, j] <- 2
    else if (cel >= quant[2, j] && cel < quant[3, j]) 
      quant.matrix[i, j] <- 3
    else if (cel >= quant[3, j] && cel < quant[4, j]) 
      quant.matrix[i, j] <- 4
    else if (cel >= quant[4, j])
      quant.matrix[i, j] <- 5
    
  }
}

# get 1st PC values
pc1 <- fit.impact$coord[,1]
# create ascending index from PC1 values (from negative to positive)
pc1.idx <- order(pc1)

# order quantile-transformed data columns by pc1 impact index
quant.matrix.order <- quant.matrix[pc1.idx]
# create matrix object from ordered quantiles by pc1 impact
matrix <- as.matrix(quant.matrix.order)

png("output/traffic-light.png", width=1900, height=1600, res=200, type="cairo")
# finally prepare pics
x <- 1:(nrow(data.noyear))
y <- 1:(ncol(data.noyear))
op <- par(mar = c(3, 3, 3, 6), oma = c(0.5,2,0,0), xpd = TRUE)
image(x, y, z = matrix, zlim = c(1,5), 
      col = c("green", "yellowgreen", "yellow", "gold", "red"), 
      axes = FALSE, xlab = "Year", ylab = "")
axis(1, at = seq(1, nrow(data.noyear), by = 1), tick = FALSE, labels = data$year, 
     cex.axis = 0.7, las = 3)
axis(2, at = seq(1, ncol(data.noyear), by = 1), tick = FALSE, 
     labels = data.names[pc1.idx],
     cex.axis = 0.8, las = 1, padj = 1) 
box()
title(main = "Traffic Light Plot", font.main = 2)
legend(x = max(x)+0.5, y = max(y), title = "Quantile colors:",
       legend = c("0 ... 0.2", "0.2 ... 0.4", "0.4 ... 0.6", "0.6 ... 0.8", "0.8 ... 1"),
       fill = c("green", "yellowgreen", "yellow", "gold", "red"),
       border =  c("green", "yellowgreen", "yellow", "gold", "red"), 
       cex = 0.7, bty = "n")
par(op)
dev.off()

# correlation plot
png("output/correlation-test-col.png", width=1900, height = 1600, res = 200, type="cairo")
cor.mat <- cor(data[,-1])
p.mat <- cor.mtest(data[,-1])$p
corrplot(cor.mat, method = "color",
         type = "upper", order = "hclust", number.cex = .7,
         addCoef.col = "black", # Add coefficient of correlation
         tl.col = "black", tl.srt = 90,
         p.mat = p.mat, sig.level = 0.05, insig = "blank", 
         diag = FALSE)
dev.off()


# prepare 3d plot by PCA dimensions 1-3
#fit.coord <- fit.score$coord[,1:3]
#fit.coord.offset <- fit.coord * 0.1
# plot3d(x = fit.coord[,1], 
#        y = fit.coord[,2], 
#        z = fit.coord[,3],
#        xlab = "Dim 1",
#        ylab = "Dim 2",
#        zlab = "Dim 3",
#        size = 10,
#        type = "p",
#        col = "#0d0e69")
# text3d(x = fit.coord[,1] + fit.coord.offset[,1], 
#        y = fit.coord[,2] + fit.coord.offset[,2], 
#        z = fit.coord[,3] + fit.coord.offset[,3], 
#        rownames(fit.coord), 
#        col = "darkgray",
#        cex = 0.7,
# )
# lines3d(x = fit.coord[,1], y = fit.coord[,2], z = fit.coord[,3], lty = "dashed")
# writeWebGL("output/3dpccoord/", width = 600, height = 800)
